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Q ■ Abstract 

\o 

r\ We present a technique which permits the calculation of two-point functions of op- 

erators containing one heavy quark and an arbitrary number of light quarks as analytic 
functions of the heavy-quark mass. It is based on the standard Jacobi linear solver 
q; used for the calculation of quark propagators. Results for the heavy- light pseudoscalar 

and vector meson masses are obtained on 16 3 x 48 lattices at (3 = 6.2 using the Wilson 
fermion action, and agree with published data. The incorporation of smeared operators 
and 0(a)- improved actions presents no problems. 



Introduction There is much interest at present in studying systems involving both heavy 
and light quarks on the lattice. The reasons for this are twofold. Firstly, the values of 
quantities such as the decay constants of the D and B mesons and the form of the Isgur- 
Wise function [|l| are currently of great interest to both the experimental and theoretical 
communities. These are, in principle, calculable on the lattice. Secondly, it is expected that 
discretisation errors are relatively large for states containing heavy quarks and such systems 
provide an arena in which to look for gains from the use of improved actions . 

Although nature has so far only revealed two heavy quarks, charm and bottom, with masses 
of approximately 1.5 and 5 GeV respectively, it is not always sufficient to perform lattice 
simulations solely at these two points. Regarding the evaluation of the Isgur-Wise function 
C,(v.v'), it is essential to know when (and if) the heavy-quark effective theory is valid. Once 
this has been determined, it is necessary to perform measurements at many different heavy- 
quark velocities to map out the form of £(i>.i>'). Since the lattice quantises the allowed 
values of momenta, this can only be done in a continuous fashion by varying the quark mass. 
In order to study lattice artifacts we need to have results at many heavy-quark masses to 
discover when and how the discretisation errors manifest themselves and to decide when to 
trust the lattice results. 

In this letter, we use the standard Wilson fermion action 



i' 



q{l- K fl)q 



where k, the quark hopping parameter, is inversely related to the bare mass. We develop 
a method which allows the hadron correlators to be evaluated as convergent power series 
in the heavy-quark hopping parameter. In what follows we denote generic quarks and their 
hopping parameters by q, q and k, and use h, h and K h for heavy quarks. 

Jacobi Algorithm To evaluate any lattice quark propagator G we need to solve the linear 
equation 

(1-kP)G = V (2) 

for G with some given source 77. The Jacobi algorithm comes from the naive Taylor expansion 

of the inverse of the fermion matrix, i.e. 

00 
G = V + J2K n $> n r}. (3) 

n=l 

In practice, it is implemented by the recursion relation 

G = V 

G n+1 = r) + K$G n (4) 



where we terminate the series at some finite n max . The minimum value of n max needed to 
achieve some desired accuracy can be chosen, for example, by monitoring the norm of the 
residue vector 

r n = T) - (1 - Kp)G n . (5) 

The algorithm is convergent for sufficiently small k, and for the heavy-quark masses of 
interest there is no problem with convergence. For our purposes it is sufficient to notice that 
at iteration n, G n is updated by a term of order K n . 

Dividing the lattice into even and odd sites and using the fact that p only connects even 
sites to odd sites it is apparent that the solution G of equation (0) satisfies 

G dd = ??odd + KpG even (6) 

We need only compute the solution on even sites which satisfies 

(1 - K 2 p 2 )G cvcn = r] cvcn + Kprjodd (7) 

Equations (H) and (|?P comprise the red-black preconditioning of equation (^) and for the 
Jacobi algorithm reduce the computational load by a factor of two. 

Hopping Parameter Expansion For simplicity, consider a mesonic operator Or = hTq. 
The quantity of interest is the timesliced correlator Cr(t) defined by 

C r (t) = ]T(0r(£,£)(^(O,O)> 

X 

= ^{Tr { l5 rG q (x,t)T l5 G{(x,t)}) (8) 

X 

where T = 7 r^7 . We assume that the computationally intensive part of this calculation 
(the evaluation of the light quark propagator G q ) has been performed. Using the Jacobi 
algorithm, the standard approach would be to solve the red-black preconditioned equation (|7|) 
for some particular value of k^ and compute Gh- Cr(t) would then be evaluated and the 
whole procedure repeated for each different value of k^- 

We introduce the hopping parameter expansion Q by inserting the Taylor series expansion 
for Gh directly into equation (H) to obtain Cr(t) as a power series in n h : 

7^ max 

Cv{t) * E 4 B) (*K (9) 

n=0 



where 



and_£9°?7 = rj. 



cP(t) = £<Tr { 75 rG,(x,t)f 75 PV (*,*)}) (10) 



The computational overhead is the calculation of the desired traces at every order of the 
expansion rather than only once for the sum of the series. The gain is that the light quark 
propagator only has to be read in to memory once, which may be significant on certain 
machines for large lattices, and results are obtained for a continuous range of Kh- 

It is obvious that the static approximation (the limit of infinite quark mass / zero Kh) is a 
simple by-product of this more general finite mass expansion. The relationship is 

Cf atic (t) = c?(t). (11) 

The numerical problem with equation (§) as it stands is that the coefficients cp (t) diverge 
as n becomes large and eventually cause an overflow on the computer. We can predict the 
growth of Cp (t) by the following simple argument. 

Consider the pseudoscalar channel (r = 75). At some critical value of the hopping parameter 
k = k c , the pseudoscalar comprising two degenerate light quarks (the pion) becomes massless. 
In this case C-p(t) no longer decays exponentially and so we expect all orders of the hopping 
parameter expansion to contribute with roughly equal weight, i.e. 

c^(t)K n c ~ constant. (12) 

Since the pseudoscalar is the lightest channel, c^(t) will have the worst divergences. If 
we can keep these under control, then all other channels will be numerically stable as well. 
Thus, although in practice we deal with light quark propagators which correspond more to 
the strange quark mass than the chiral limit, if we rewrite equation (|9|) as 

^max 

c r (t)~ Et'WH" (is) 



n=0 



»/ 



with a ~ k c 1 then we expect the new expansion coefficients of (t) to be roughly constant. 

Simulation Details Numerical studies were performed on 10 lattices of size 16 3 x 48, 
generated in the quenched approximation at (3 = 6.2. Since the major purpose of this letter 
is to demonstrate the validity of the technique and to reproduce known results, only two 
meson operators were calculated: 



Pseudoscalar : 


r = 


75 


Vector : 


r = 


3 
£7, 



from which the masses and hyperfine splitting can be calculated and compared to reference 
Q. The light-quark hopping parameter was chosen as k = 0.1510, corresponding roughly 
to the strange quark mass ||, and these light propagators were calculated using an over- 
relaxed minimal residual algorithm. To span the range of heavy- quark masses in |4|] (0.125 < 
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Figure 1: Convergence of Vector (Cy) and Pseudoscalar (Cp) Meson Correlators 



Kh < 0.145), the hopping parameter expansion was computed to 200 th order. This number 
was chosen by explicitly calculating quark propagators at Kh = 0.145 using the Jacobi 
algorithm and requiring that the residue vector of equation (^) satisfy |r nma J < 10~ 7 . The 
following test was performed on a single configuration. Firstly, the standard propagator code 
(used to evaluate the light propagator) was used to compute the heavy-quark propagators 
at Kh = 0.145 and Kh = 0.125, and the values of Cr(t) calculated directly. Secondly, the 
coefficients Cp (t) were computed and the series summed off-line at the same two values of 
Kh- The results are illustrated in figure [I], where the solid line is the hopping parameter 
expansion summed to a given order and the dashed line is the value of Cr(t). The results of 
these two independent calculations agreed, and it was deduced that the hopping parameter 
expansion had converged (to the accuracy desired) within the required range. As expected, 
convergence is most rapid at small values of Kh and t. The other notable features are that 
the convergence of neither the vector nor the pseudoscalar meson correlators is monotonic, 
and that the full 200 orders are necessary at the lightest mass and furthest timeslice. For 
heavy quarks near the charm mass (Kh — 0.135 0) around 100 orders are sufficient. 

The value of a was varied, and the scaled coefficients cf 1 (£) found to be numerically stable 
for a = 6.0. The bare coefficients c^(t) actually grow like (6.3) n which, since the critical 



hopping parameter k c ~ 0.1533 ||, agrees with our prediction c^(t) ~ (6.5) n . 

For these initial studies, point sources at the origin (an even site) and point sinks were used 
for both the heavy and the light quarks. In this case, the hopping parameter expansion is 
an implementation of the red-black preconditioned Jacobi algorithm. Although our tests, 
performed on a single configuration, were sufficient to prove that the algorithm works we 
decided to use smeared operators during the full simulation in order to obtain better signals 
for the relevant particles. Point sources were still used for the heavy quarks so that the 
red-black factorisation could be exploited. The light-quark operators were smeared at the 
source and the sink using a gauge covariant smearing function with an RMS smearing radius 
of about four lattice spacings ||. The incorporation of non-local heavy-quark sources and 
sinks is discussed later. 

Results Since in principle the simulation provides Cr{t) as an analytic function of Kh 
there are several analysis methods available. Previously, when the expansion was truncated 
well before convergence, the behaviour of the Fourier transformed correlation functions was 
studied and the particle masses identified as poles in momentum space @. Our approach is 
to sum the series at many values of Kh and fit the asymptotic exponential decay in time. 

After the use of smeared operators, plateau regions in Cr(t) were readily identified for the 
whole range of Kh, although the limited statistics meant that the signal deteriorated rapidly 
especially at very large heavy-quark masses (here we approach the static quark limit where 
this effect is well known ||). An uncorrelated fit to the vector and pseudoscalar meson chan- 
nels was done from time slices 5 to 10 and the corresponding masses extracted. To estimate 
the errors, the whole procedure was bootstrapped and the probability distributions of the 
masses and mass difference obtained. The quoted errors are 68% confidence level regions 
and are therefore not necessarily symmetric. Unfortunately, lack of statistics prevented a 
proper correlated fit. 

The most interesting quantity is my — m 2 P , which becomes constant in the heavy-quark 
effective theory with corrections of order (my + rrip)' 1 . Experimentally, this quantity is 
remarkably constant all the way from the n, p system through the K and D to the B, 
with a measured value of around 0.55 GeV 2 . We plot this mass difference against m 2 P in 
figure H] together with the data from ||. Our analysis was actually done at 17 evenly spaced 
values of Kh, and we have drawn a continuous line representing the best fits and dashed 
lines representing the envelope of the error bars. It is clear that, within errors, there is 
no discrepancy between the two data sets. The results show the well known feature that 
the hyperfine splitting for heavy-quark mesons is far too small for the Wilson action (at 
(3 = 6.2, a -1 = 2.73 GeV ||) and decreases with increasing heavy-quark mass, contrary to 
experiment. This is observed both in heavy-light and heavy-heavy systems. 
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Figure 2: Vector-Pseudoscalar Hyperfine Splitting 

Non-Local Operators and Improved Actions It is not practical to smear the heavy- 
quark propagator at the sink since this would have to be done at each order of the hopping 
parameter expansion, taking a prohibitively long time for even the simplest of smearing 
techniques. If smearing is required at the sink, it should be performed on the light quark 
fields as this only needs to be done once. The only complication arises if smearing is required 
at the source because the only light-quark propagators available originate from point sources. 
In this case, it is necessary to smear the heavy quark at the source. There is no problem in 
doing this, the only disadvantage being that the red-black factorisation cannot be exploited. 

It has been suggested that the unphysically small value obtained for the hyperfine splitting 
in heavy-quark systems is due to large order a effects present in the Wilson action. Using 
an 0(a)-improved fermion action would be expected to improve the situation, and this has 
been demonstrated for the "clover" action by studying the charmonium system JFL M. The 
important point is that k still only enters the fermion action as an overall prefactor so it is 
straightforward to use the hopping parameter expansion. The other element of the clover 
action is that the quark propagator needs to be rotated at source and sink to complete the 
0(a) improvement 0, so there is no advantage to red-black factorisation. 



Conclusions We have shown that the hopping parameter expansion is a viable method for 
the study of heavy-light systems on the lattice. It provides results as an analytic function of 
the heavy-quark mass provided that this mass is larger than some lower bound. This lower 
bound can be reduced by going to higher orders of the expansion, but for studies at (3 — 6.2 



on a 16 3 x 48 lattice only 200 orders are needed for Kh < 0.145. This is lighter than charm 
for which Kh — 0.135 and where around 100 orders are sufficient. The hopping parameter 
expansion, when used with an 0(a)-imp roved action, provides an attractive method for map- 
ping out the behaviour of meson decay constants and the Isgur- Wise function continuously 
at and above the charm quark mass. 
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